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Turbulence is a remarkable subject in physics. The underlying equations, which 
are in their simplest formulation the Euler equations, were published 250 years ago 
Il20l . Yet a theoretical grasp of the phenomenology emerging from these equations had 
not been achieved before the mid-twentieth century, when Heisenberg 11251 and Kol- 
mogorov [33 1 obtain first analytical results. Eventually, it took the capabilities of mod- 
ern supercomputers to obtain a full appreciation of the complexity that is inherent to the 
Euler equations. Astrophysics is now at the very frontier of numerical turbulence mod- 
elling. Among the additional ingredients for making turbulence in astrophysics even 
more complex are supersonic flow, self-gravity, magnetic fields and radiation transport. 
In contrast, terrestrial turbulence is mostly incompressible or only weakly compressive. 
External gravity is, of course, an issue in the computation of atmospheric processes on 
Earth. Self-gravity, however, is only encountered on large, astrophysical scales. The 
dynamics of turbulent plasma has met vivid attention in research related to nuclear 
fusion reactors but, otherwise, is not encountered under terrestrial conditions. 

In this Chapter, I give an overview of the various approaches toward the numerical 
modelling of turbulence, particularly, in the interstellar medium (ISM). The discussion 
is placed in a physical context, i. e. computational problems are motivate from basic 
physical considerations. Presenting selected examples for solutions to these problems, 
I introduce the basic ideas of the most commonly used numerical methods. For de- 
tailed methodological accounts, the reader is invited to follow the references. Some 
important results and astrophysical implications are briefly outlined. Since turbulence, 
in a strict sense, is genuinely three-dimensional, I almost exclusively consider three- 
dimensional simulations. Furthermore, turbulence is a multi-scale phenomenon and, 
in order to capture its properties correctly, sufficient numerical resolution is essential. 
This is why newer higher-resolution computations are generally preferred as exam- 
ples (as regards self-gravitating turbulence, however, high-resolution simulations are 
appearing just now and are not included here). 

To begin with, I briefly discuss the equations which are numerically solved in as- 
trophysics and related theoretical aspects. The following Sections deal with supersonic 
turbulence, self-gravitating turbulence and magnetohydrodynamic turbulence. Natu- 
rally, there are intersections. Representative examples for numerical simulations are 



1 



placed depending on their main objectives. Finally, I give a resume of what has been 
achieved and which challenges are likely to be met. 



1 Fundamentals 

In the following, we consider self-gravitating inviscid fluid subject to the action of 
external force fields and, possibly, magnetic fields. The evolution of the mass density 
p, the velocity v and the specific energy e of the fluid is given by the compressible Euler 
equations: 

gp = -pv», (1) 
p^-v = -vp+p(f + g), (2) 

p—e + VvP = r-A + pv(f + g), (3) 
where the Lagrangian time derivative is defined by 

— = — + v ■ V. (4) 

Dt at 

The total energy per unit mass is given by 

e -V + r-V< (5) 

2 (y-l)p 

where y is the adiabatic exponent and the pressure P is related to the mass density and 
the temperature T via the ideal gas law: 

P = P -^. (6) 

The constants kg, fi and ma denote, respectively, the Boltzmann constant, the mean 
molecular weight and the mass of the hydrogen atom. The energy budget can be altered 
by heating (T) and cooling (A), by non-gravitational forces (pf) as well as gravity 
ipg). Generally, heating and cooling processes are important for the dynamics of the 
interstellar medium (Sections[3]&0. 

Non-gravitational forces supplying energy to the flow can be mechanical or mag- 
netic. An example of a mechanical force would be the random driving force that is 
commonly used in turbulence simulations. This is an external force, i. e., it is indepen- 
dent of the dynamics of the fluid (Section|2]l. Quite the opposite holds for conducting 
fluids in the presence of magnetic fields. In the case of ideal magnetohydrodynamics, 
the fluid is dragged by the Lorentz force pf — J x B, where the current / = V x B. 
The magnetic field B, in turn, depends on the flow via Faraday's law, 

dB 

— = V x (v x B). (7) 

at 
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The interaction between turbulent flow and magnetic fields alters properties of turbu- 
lence such as the scaling behaviour (Section|Ul. 

The gravitational acceleration g = -V ■ <p arises, in part, from the self-gravity of 
the fluid and, depending on the boundary conditions, the gravitational field of exter- 
nal sources. In the case of periodic boundary conditions, gravity is solely produced 
by density fluctuations with respect to the global mean (Section [3j. In this case, the 
gravitational potential is determined by the Poisson equation 

V 2 = 4ttG( P - <p», (8) 

where G is Newton's constant. The mean mass density (p) = po is a constant because 
of mass conservation. 

Letting the curl operator Vx act upon equation (0, an evolutionary equation for the 
vorticity of the flow, to = V x u, is obtained: 

D 1 

— fc> = S • GJ - + ^rVp X VP + V X /. (9) 

Dt p 2 r 

The rate-of-strain tensor S is the symmetrised velocity gradient and the divergence 
d — V • v. Apart from the stirring of fluid due to rotational force components, vorticity 
can be generated by the baroclinic term Vp x VP, if the gradient of the mass density 
and the pressure gradient are not aligned. With the ideal gas equation [6] it follows 
that baroclinic vorticity generation occurs in non-isothermal gas. The term S ■ a> - dco 
accounts for the stretching of vortices caused by strain in addition to the expansion or 
contraction of vortices due to the compressibility of the fluid. 
The time evolution of the divergence d is given by 

p 2 

This equation also follows from the conservation law <J2J by applying the operator V- 
and substituting the Poisson equation for the gravitational potential ([8j. The norm of 
the rate of strain is defined by |5| = (25 jjS ij) 1 ^ 2 - From the first term, one can see that 
vorticity contributes to positive divergence, while strain is associated with negative 
divergence, i. e. converging flow. In the incompressible case, the first term on the right 
hand side of equation ( TlOb is cancelled by the Laplacian of the pressure divided by 
the density. Whereas vorticity cannot be generated by a gradient field such as gravity, 
an increasing rate of convergence will be caused by gravity, if the term 4nG(p - po) 
dominates. In this case, gravitational collapse of the gas ensues. 

The numerical simulation of astrophysical turbulence has to account for all fluid 
dynamical processes which become manifest in the evolutionary equations of the vor- 
ticity © and the divergence (ITOb . In many terrestrial applications, on the other hand, 
incompressible turbulence without gravity and magnetic fields is considered. Then 



1 9 

-V 2 P 
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-47TG(p-p ) + V-/. 



(10) 
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equations (0 and ( TTOb reduce to 



D 




S ■ co + V x /, 



(11) 



0. 



(12) 



The only possibility of generating vorticity in this case is to stir the fluid by solenoidal 
(rotational) forces. Non-linear turbulent interactions stretch and fold the large-scale 
eddies produced by stirring into thin vortex filaments. This is the very essence of tur- 
bulent fluid dynamics. Here, the question arises whether the generation of turbulence 
in the ISM resembles the scheme of stirring, stretching and folding at all. Interest- 
ingly, this presumption underlies the application of solenoidal driving forces in most 
numerical simulations. 

Incompressible fluid dynamics serves as an important reference case for which scal- 
ing properties of isotropic turbulence in the ensemble average can be derived analyti- 
cally. Kolmogorov IT331 showed that the root mean square velocity fluctuations u'({) at 
a certain length scale I obey the so called 2/3-law, 



provided that £ is small compared to the integral length scales L at which energy is 
supplied to the flow. The 2/3-law implies that the rate of kinetic energy dissipation 
e is asymptotically constant, as e ~ v'(€f ji =s const, in the limit £/L «: 1. This 
result is known as the law of positive energy dissipation 11221 . Both the 2/3 law and 
the law of positive energy dissipation appear to be robust properties of turbulence as 
long as the flow becomes isotropic and nearly incompressible towards small length 
scales. This is not always the case. Magnetic fields, for instance, introduce small-scale 
anisotropy. Moreover, there is no consensus yet to what extent the scaling properties 
of incompressible turbulence carry over to supersonic flow. 

2 Supersonic turbulence 

Turbulence is called supersonic, if the root mean square (RMS) Mach number A1 rms = 
((u/cs) 2 ) 1 / 2 exceeds unity. Kinetic energy dissipation in supersonic turbulence caused 
by shock fronts propagating through the fluid dominates viscous dissipation. Since the 
resulting dissipative heating rises the speed of sound to a value comparable to the typi- 
cal flow velocity within a dynamical time, supersonic turbulence cannot be maintained, 
if the excess heat is not efficiently removed from the fluid. As an idealisation, heating 
is exactly balanced by cooling such that the temperature remains constant. In this case, 
continuously driven turbulent flow settles into a steady state in which the RMS Mach 
number is asymptotically constant. In isothermal gas, solenoidal forcing is required to 
produce turbulence because the density gradient is aligned with the pressure gradient 
and, consequently, Vp x VP = 0. In this case, it can be seen from equation (0 that 
non-zero vorticity either must be imposed as initial condition or it is generated by the 
curl of the force field. As a consequence of the constant speed of sound, the properties 
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of driven supersonic turbulence in isothermal gas depend on a single parameter only, 
namely, the RMS Mach number (provided that the forcing is purely solenoidal l50lD . 

A a lot of effort has gone into the numerical simulation of supersonic flow by means 
of finite-volume schemes of higher order, in particular, the piecewise parabolic method 
(PPM) 0171 . The applicability of PPM to turbulence was tested in simulations of adi- 
abatic turbulence [53 1 and nearly isothermal transonic turbulence. An example for the 
latter is the driven turbulence simulation by Porter & Pouquet ll47l . They utilised an 
explicit cooling function and, thereby, maintained a root mean square Mach number 
~ 1.0. In ||49l , the equation of state of degenerate electron gas was applied in simula- 
tions of forced transonic turbulence with PPM. Degeneracy entails an extremely large 
heat capacity. Since the isothermal case corresponds to the limit of infinite heat capac- 
ity, the gas in these simulations is effectively isothermal. One of the conclusions was 
that the numerical dissipation of PPM acts as an implicit subgrid scale model in the 
weakly compressible regime (at small scales). 

A new approach is the application of adaptive mesh refinement (AMR) in turbu- 
lence simulations [36 1. AMR is based on the finite-volume approach with a hierarchy 
of grid patches of different resolution. The basic idea is to represent relatively smooth 
flow regions on coarse grids only, whereas steep gradients of the velocity, density or any 
other field are computed with higher accuracy by dynamically inserting grids of higher 
resolution j6]|7]. In astrophysics, AMR has been successfully used for the treatment of 
strong shocks or gravitational collapse in a variety of problems [44 1 . On the other hand, 
the simulation of developed turbulence with AMR has been regarded as infeasible, be- 
cause, in the picture of the turbulent cascade of eddies, small-scale features of the flow 
would be space-filling [ 22 1 . For a particular flow realization, however, turbulent eddies 
occupy an ever decreasing fraction of the flow domain toward smaller length scales at 
any given instant of time. This follows from the intermittency of turbulence 11221 . 

Picking up the intermittent picture of turbulence, Kritsuk et al. proposed that AMR 
would offer a computational advantage compared to static grids of fixed resolution 11361 . 
This can be motivated from the Landau estimate of the number of degrees of freedom 

Eg, 

N ~ Re 9/4 , (14) 

where Re is the Reynolds number of the flow. Present day supercomputers manage N ~ 
10 10 , which allows for Re ~ 10 4 . Note that fully developed turbulence is known to set 
in at Reynolds numbers of a few thousand. On the other hand, invoking intermittency 
models of turbulence such as the yS-model 11221 . we have 

N „ Re^+D, (15) 

where D is interpreted as the fractal dimension of dissipative structures. Assuming a 
value of D « 2 as in the Burger model of supersonic turbulence, it follows that N ~ 10 8 . 
This figure is smaller than the Landau estimate (corresponding to a fixed-resolution 
grid) by two orders of magnitude. Of course, there is substantial computational over- 
head with AMR in comparison to static grids. Apart form that, more sophisticated in- 
termittency models and numerical studies suggest that D > 2 in supersonic turbulence 
IfTOllTTI . Nevertheless, AMR is expected to pay off, if several levels of refinement are 
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Figure 1 : Projected mass density in an AMR simulation of supersonic turbulence with 
effective resolution 2048 3 by Kritsuk et al. 2007, using Enzo lf37l . Bright regions 
contain gas of high density, while under-dense gas appears dark. 

used. In particular, this applies to scenarios including gravitational collapse (see the 
next Section). 

In |[37ll , the effective resolution (i. e. the resolution corresponding to the most 
refined grids) of an AMR simulation of supersonic turbulence was raised to the yet 
unprecedented value of 2048 3 . Nearly isothermal gas was maintained by setting the 
adiabatic exponent y to a value that differs only by a small fraction from unity. This 
means that the internal energy E- mi = P/(y— l)p is artificially increased by a huge factor 
such that E mt » |pi> 2 and, hence, it takes many dynamic time scales to heat the gas 
significantly by kinetic energy dissipation. In combination with a 1024 3 static grid sim- 
ulation, a great wealth of results was obtained from the analysis of the numerical data. 
Figure Q] shows the projected mass density in logarithmic scaling for one snapshot of 
the AMR simulation. One can see voids in between high-density regions which display 
intricate turbulent structure. This is a tell-tale sign of the pronounced intermittency of 
supersonic turbulence. 

Further indications of intermittent properties come from the probability density 
function of the time-averaged mass density, which follows very closely a log-normal 
distribution |45|, and the scaling exponents of the velocity structure functions. The 
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definition of structure functions is based on spatial correlations of the velocity field: 

s p (e) = (\v(x + ()-v(x)\y. (16) 

The square brackets denote the average over the whole domain of the flow. Structure 
functions probe the correlations of the velocity field at different spatial positions de- 
pending on the separation I. It is both a theoretical prediction and an experimentally 
well established fact that the structure functions of isotropic turbulence obey power 
laws 

S p (£)k&. (17) 

The exponents ( p are called the scaling exponents. For incompressible turbulence, Kol- 
mogorov found ( p = p/3. In the case of the second-order structure functions (p-2) this 
result corresponds to the 2/3-law for the velocity fluctuations mentioned in SectionQ] 
Calculating the scaling exponents form their numerical data, Kritsuk et al. were able 
to demonstrate significant deviations form the Kolmogorov relation in accordance with 
predictions of the intermittency model proposed by Boldyrev et al. ATI [37ll , which 
is a generalisation of the She-Leveque model for incompressible turbulence 11481 . In 
particular, it was found that £2 ~ 1-0, whereas (2 = 2/3 in the Kolmogorov theory. 
Intriguingly, the analysis also revealed that the scaling exponents applicable to incom- 
pressible turbulence are obtained, if the statistics is computed for the density-weighted 
variable p l ^ 3 v in place of i>, a result that is not fully understood yet but might bear 
important implications on the nature of supersonic turbulence. 



3 Self-gravitating turbulence 

According to the linear stability analysis of the Euler equations by Jeans [28], a per- 
turbation in isothermal gas of temperature T and uniform density po becomes unstable 
against gravitational collapse, if its size exceeds the Jeans length 

Ai = c *\[t^-> (18) 

where c s oc T 1 ^ 2 is the isothermal speed of sound. Beginning with Chandrasekhar's 
proposition to substitute c 2 by 

4 =c 2 s + \(» 2 ) (19) 

in the above expression for the Jeans length, several attempts have been made to ex- 
tend the Jeans stability analysis to turbulent gas of RMS velocity (v 2 ) 1 ^ 2 lfl6l fl2l FPU . 
Equation (T% implies that the effective gas pressure is given by the sum of the thermal 
pressure P and the turbulence pressure P tur b = l-Ekin, where is the mean kinetic 
energy density. 

However, even the most advanced analysis put forward so far lfT3l . suffers from 
severe constraints. A perturbation analysis can be carried out for a statistically sta- 
tionary equilibrium state only. Thus, it has to be assumed that the free fall time scale 
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Iff ~ (Gpo) is much greater than dynamical time scale of turbulent flow. From the 
divergence equation ( [Tol l, one can see that this assumption implies that the self-gravity 
term remains small compared the other terms at all times. Allowing for gravitational 
collapse, however, self-gravity eventually dominates and 

^■d ~ -Gp. (20) 

A non-perturbative theory of the regime in which the free-fall time scale and the dy- 
namical time scale of turbulence are comparable was suggested by Biglari and Dia- 
mond |8 |. Combining scaling relations from the yS-model of turbulence ll22l with the 
assumption of energy equipartition between gravity and turbulence at all scales, they 
derived an intermittent hierarchical cloud model of self-gravitating turbulence. 

At present, however, there is no general theory that would fully encompass the 
highly non-linear and non-stationary interplay between gravity and turbulence in the 
interstellar medium. Numerical simulations, on the other hand, can aid to the under- 
standing self-gravitating turbulence, although one has to resort to artificial mechanisms 
of producing turbulence or imposing more or less arbitrary initial conditions. One tech- 
nique is to apply a driving force as mentioned in the previous Section. Alternatively, 
some random initial velocity field might be assumed as initial condition. This results in 
decaying turbulence. Yet another option is to consider gas in an unstable state and apply 
small perturbations. The most prominent example is the thermal instability ll35l[T9l . 

Numerical simulations of self-gravitating driven or decaying supersonic turbulence 
were initially performed with smoothed particle hydrodynamics (SPH) because of its 
ability of tackle variations of the mass density over many orders of magnitude BP . 
SPH makes use of the Lagrangian framework of fluid mechanics and represents fluid 
parcels by particles. Examples for the application of SPH to self-gravitating turbulence 
are the simulations by Klessen et al. l|29l l30l |3D . In these simulations, the influence 
of self-gravity on probability density functions 11291 and the fragmentation properties 
of the gas was investigated both for driven and for decaying supersonic turbulence 
11301 |3~T1 . Sink particles were introduced to capture collapsing regions beyond a cer- 
tain density threshold in order to prevent the mass density from growing indefinitely. 
Figure |2]illustrates that sink particles are mostly formed within filamentary structures 
(which are also seen in cosmological simulations). Varying the Mach number, it was 
found that higher turbulence intensity slows down the rate at which sink particles are 
produced considerably (Figure [3). The major conclusion drawn by Klessen et al. was 
that supersonic turbulence inhibits gravitational collapse globally, as one would sus- 
pect from equation d!9t . The strong gas compression caused by shocks, on the other 
hand, can trigger gravitational collapse locally and initiate the formation of stars. 

Jappsen et al. Il27l studied the mass spectrum produced by turbulent fragmentation 
in non-isothermal gas. They defined the state variables by piecewise polytropic rela- 
tions, i. e. P oc pW for a certain range of densities. If y e ff < 1, the gas cools with 
increasing density. This effect enhances the compressibility of the gas and therefore 
supports gravitational collapse (negative Vp • VT contribution adds to the gravity term 
in the divergence equationlTOli . At some threshold density, however, the cooling regime 
ceases and y e g becomes greater than unity. From various simulations, it was found that 
the transition from the cooling regime to isothermal or nearly isothermal gas sets a 
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Figure 2: Filamentary structure with sink particles in an SPH simulation of self- 
gravitating turbulence (courtesy of R. Klessen). In this plot, over-dense regions appear 
bright. 

characteristic mass scale of turbulent fragmentation. This mass scale was interpreted 
to correspond to the peak of the observed initial mass function. Vazquez-Semadeni 
et al. ll55ll included explicit heating and cooling in SPH simulations of colliding gas 
streams. They made use of a model for the cooling function A in the energy equa- 
tion (0 that was proposed by Koyama & Inutsuka [34|. Due to the thermal instability 
(compression causes the gas to cool), the gas at the collision interface undergoes grav- 
itational collapse and becomes increasingly turbulent. This is interpreted as possible 
mechanism of molecular cloud formation. Remarkably, approximate equipartition of 
gravitational and kinetic energy was found, although clearly no state of virial equilib- 
rium was approached by the collapsing gas. Further applications of SPH addressing 
problems such as the formation of brown dwarfs, binary star systems and stellar clus- 
ters via turbulent fragmentation as well as the origin of the initial mass function were 
presented by Bate et al. OIHEI an d Bonnel et al. lfT4ll . 

Despite the numerous important contributions to the understanding of the role of 
turbulence to star formation, the adequacy of treating turbulence with SPH has been 
a matter of debate. On the one hand, numerical studies of self-gravitating systems 
indicate basic agreement between AMR and SPH [18|. However, comparisons with 
grid-base codes suggest that SPH dissipates small-scale velocity fluctuations signif- 
icantly stronger ll46l [32l . In particular, it was demonstrated that the steeper power 
spectra obtained from SPH simulations are accompanied by mass distributions which 
are not consistent with observations. Apart from that, it appears to be rather difficult 
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Figure 3: Total mass fraction captured by sink particles as a function of time in 
SPH simulations with different Mach numbers and varying driving scale of turbulence 
(Klessen et al., 2000 flH). 

to accommodate magnetohydrodynamics in the SPH formalism. Although there are 
ongoing attempts to overcome this shortcoming |40|, MHD has been treated with great 
success using finite-volume discretisation. 

4 Magnetohydrodynamic turbulence 

In purely hydrodynamical turbulence, there is no preferred spatial direction for small- 
scale velocity fluctuation. The randomisation of the flow due to non-linear turbulent 
interactions produces statistical isotropy towards small scales which can be interpreted 
as a symmetry of the ensemble average ll22l . This symmetry is broken by a magnetic 
field B in turbulent conducting fluid, if the ratio of the thermal to the magnetic pressure, 

is comparable to unity or less. If the curl of the Lorentz force, V x (/ x B), in the 
vorticity equation (O is sufficiently strong, then the vorticity will mainly grow in the 
direction of B. In developed magnetohydrodynamic (MHD) turbulence, this effect 
causes eddies to shrink perpendicular to the field lines. As a further implication, the 
anisotropy of MHD turbulence increases towards smaller length scales |9). 

The numerical simulation of MHD turbulence is rather challenging. Godunov- 
based schemes of higher order such as PPM become very complex upon including 
MHD. For this reason, simpler schemes have been adopted for MHD. Examples are 
the Zeus code l52l . the Stagger code l43l . and the Ramses code which also features 
AMR (541 [231 . As an illustration, Figure [4] shows a simulation of the gravitational 
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Figure 4: Simulation of a collapsing gas cloud of 1000 solar masses with Ramses 
(courtesy of P. Hennebelle). The magnetic energy of the clouds is about the same as its 
kinetic energy. 

collapse of a magnetised cloud with up to 6 levels of refinement. One of the major 
problems that needs to be addressed when solving the MHD equations numerically is 
to keep the magnetic field divergence-free, i. e. V ■ B = 0. In Ramses and other codes, 
the constrained transport algorithm is employed to ensure this constraint ll2D . 

A comparison of the properties of driven hydrodynamic (HD) and MHD turbulence 
in high-resolution simulations was presented by Padoan et al. [46 1. Whereas the energy 
spectrum functions, irrespective of magnetic fields, were found to follow closely the 
power law ]C 2 in the inertial subrange, the fragmentation properties of the gas appear 
to be markedly different for HD and MHD turbulence, respectively. In plots of the 
projected density fields for both cases (Figures[5]and[6]l, one can see that the gas is more 
concentrated in pronounced filaments under the action of magnetic fields. Utilising a 
clump-find algorithm, it was demonstrated that a significantly steeper mass spectrum of 
dense cores (defined by a density threshold based on the Bonnor-Ebert mass) resulted 
from purely hydrodynamical turbulence. Furthermore, they were able to reproduce 
the Chabrier initial mass function |[T5l in the MHD case, although this result must be 
considered with care given the ambiguity of calculating mass spectra. 

Glover and Mac Low E4l performed simulations of self-gravitating MHD turbu- 
lence including various thermal and chemical processes. They were able to show that 
the fast formation of molecular hydrogen within a few Myrs, which is typically ob- 
served in molecular clouds, can be explained as an effect of turbulence. Hennebelle 
& Audit [26 1 investigated thermally bistable MHD turbulence and found mass spectra 
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Figure 5: Projected mass density in a simulation of hydrodynamic turbulence with the 
Stagger code (Padoan et al., 2007 |46|). For comparison with MHD turbulence, see 
Figure [6] 



similar to those inferred from CO observations of molecular clouds. Although first 
computed in only 2 dimensions, the setup was generalised to include self-gravity and 
computed in 3 dimensions as well. The ansatz by Hennebelle et al. is different from 
simulations of isothermal MHD turbulence as it does not rely on a driving force or an 
initial turbulent velocity field. 

5 Perspectives 

Up to now, most approaches to the numerical simulation of astrophysical turbulence 
have focused on particular aspects, be it a thorough understanding of supersonic tur- 
bulence, the study of the dynamics of self-gravitating turbulent gas or an emphasis on 
the effects of magnetic fields. While impressive progress has been made in the area 
of isotropic turbulence, both in the HD and the MHD case, more complex scenarios 
remain elusive. The treatment of thermal and chemical processes is still fairly approx- 
imate. Our understanding of self-gravitating turbulence remains poor, despite many 
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Figure 6: Projected mass density in a simulation of magnetohydrodynamic turbulence 
with the Stagger code (Padoan et al., 2007 [46 1). For comparison with purely hydrody- 
namic turbulence, see Figure [5] 



ongoing efforts. Radiation HD and MHD in combination with adaptive methods, are 
just beginning. 

There is one aspect that should be highlighted whenever we are talking about the 
simulation of turbulence in the interstellar medium. The huge range of length scales 
from the galactic disk down to proto-stellar cores does not allow for coverage by a 
single numerical simulation. In general, there are two numerical cutoffs, one toward 
larger scales, the other toward smaller scales. The former poses the question of initial 
as well as boundary conditions and whether simple models such as periodic boxes and 
random forcing can account for those in a self-consistent fashion. The small-scale 
cutoff would, in general, necessitate closures. These are considered as essential in 
atmospheric sciences but have met very little attention in astrophysics. 

As regards the large scales, it certainly will not be feasible in the near future to run 
a simulation of a disk galaxy or even a piece of a galactic disks including all processes 
of turbulence production, the different phases of the interstellar medium, the formation 
of molecular clouds and the collapse of gas giving rise to the birth of stars. If energy is 
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Figure 7: Projected mass density in 128 3 simulations of turbulence driven by a 
solenoidal and a compressive force field, respectively (courtesy of C. Federrath). 



injected via radom forcing in a simulation of turbulence over some subrange of scales in 
the ISM, one faces the question of what the influence of the forcing might be. That the 
outcome can vary depending on the applied forcing is illustrated by a simple numerical 
experiment. Figure [7] shows projected mass densities obtained from 128 3 simulations, 
in which purely solenoidal (left) and compressive (right) forcing was applied. Although 
there is no fully developed turbulence at such low resolution, the plots suggest that 
the flow morphology is genuinely different, with markedly higher density contrasts 
and more pronounced intermittency in the case of compressive forcing. Indeed, this 
conclusion is confirmed by high-resolution simulations BP . 

For the high-fidelity treatment of small scale effects, AMR emerges as the most 
promising tool, simply following the maxim to go to very high levels of refinement 
wherever the critical events take place. This was most impressively executed by Abel et 
al. fl] in AMR simulations of the formation of the first stars in the Universe. However, 
circumstances are not that favourable when it comes to galactic star formation and 
it remains to be seen whether AMR by itself will hold its promises. An alternative 
approach was suggested by Niemeyer et al. l42l . The basic idea is to combine the 
techniques of AMR, which have been brought to great success in astrophysics, and 
subgrid scale modelling, which is commonly used in engineering. Furthermore, the 
incorporation of additional physics such as thermal processes, complex chemistry as 
well as radiation transport in AMR simulations are essential for realistic models of 
turbulence in the ISM. 
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